{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RMP3 与 RMP4 能量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> 创建时间：2019-11-01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这一节我们讨论 RMP3 与 RMP4 的能量计算。\n",
    "\n",
    "读过 Szabo 第六章的人相信对 MPn 方法有所了解。不过我们这里不关注 MPn 方法的公式推导，并且将眼光局限于 Restricted 方法。事实上写这篇文档时，尽管曾经学习过一般的 MPn 与 CCPT 的推导方式，但并没有尝试推导过 Restricted 情况下的推导，仅仅是将书上出现的公式程序化而已。\n",
    "\n",
    "这一篇文档的主要参考是 Helgaker et al. 2013 教材 [^Helgaker-Jorgensen.Wiley.2013] 的 section 14.4。关于 Spin-Orbital MP3 与 RMP3 的另一种实现方式，可以参考以 Szabo [^Szabo-Ostlund.Dover.1996] 公式为蓝本的 Psi4NumPy 简要代码 [^psi4numpy-mp3]。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from pyscf import scf, gto\n",
    "\n",
    "from functools import partial\n",
    "np.einsum = partial(np.einsum, optimize=[\"greedy\", 1024 ** 3 * 2 / 8])\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import re\n",
    "\n",
    "np.set_printoptions(5, linewidth=150, suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 分子体系与标准结果"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们所使用的分子是非对称的双氧水分子，基组为 6-31G。比较常用但名称不太不常见的变量有\n",
    "\n",
    "- `so`, `sv` 表示占据轨道、非占轨道的分割 (split)\n",
    "\n",
    "- `C`, `e` 分别表示分子轨道系数 $C_{\\mu p}$ 与轨道能 $\\varepsilon_p$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol = gto.Mole()\n",
    "mol.atom = \"\"\"\n",
    "O  0.0  0.0  0.0\n",
    "O  0.0  0.0  1.5\n",
    "H  1.0  0.0  0.0\n",
    "H  0.0  0.7  1.0\n",
    "\"\"\"\n",
    "mol.basis = \"6-31G\"\n",
    "mol.verbose = 0\n",
    "mol.build()\n",
    "\n",
    "natm = mol.natm\n",
    "nao = nmo = mol.nao\n",
    "nocc = mol.nelec[0]\n",
    "nvir = nmo - nocc\n",
    "so, sv = slice(0, nocc), slice(nocc, nmo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-150.5850337808384"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scf_eng = scf.RHF(mol)\n",
    "scf_eng.conv_tol = 1e-12\n",
    "scf_eng.max_cycle = 128\n",
    "scf_eng.kernel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "C, e = scf_eng.mo_coeff, scf_eng.mo_energy\n",
    "Co, Cv = C[:, so], C[:, sv]\n",
    "eo, ev = e[so], e[sv]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们的参考值来自于 Gaussian 的计算 (输入卡 {download}`H2O2-MP4.gjf`、输出文件 {download}`H2O2-MP4.out`)。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"H2O2-MP4.out\", \"r\") as f:\n",
    "    gaussian_output = f.readlines()\n",
    "\n",
    "def gaussian_line(string):\n",
    "    for line in gaussian_output:\n",
    "        if string in line:\n",
    "            return line[:-1].replace(\"D+\", \"E+\").replace(\"D-\", \"E-\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中一些结论有："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MP2 Corr:    -0.2690117593\n",
      "MP3 Corr:     0.0074517392\n",
      "MP4 S   :    -0.0043821604\n",
      "MP4 D   :    -0.0088190269\n",
      "MP4 T   :    -0.0067126030\n",
      "MP4 Q   :     0.0011355775\n"
     ]
    }
   ],
   "source": [
    "ref_mp2_corr = float(gaussian_line(\"E2 =\").split()[2])\n",
    "ref_mp3_corr = float(gaussian_line(\"E3=\").split()[1])\n",
    "ref_mp4_S = float(gaussian_line(\"MP4(S)=\").split()[1])\n",
    "ref_mp4_D = float(gaussian_line(\"MP4(D)=\").split()[1])\n",
    "ref_mp4_Q = float(gaussian_line(\"MP4(R+Q)=\").split()[1])\n",
    "ref_mp4_T = float(gaussian_line(\"E4(SDTQ)=\").split()[1]) - float(gaussian_line(\"E4(SDQ)=\").split()[1])\n",
    "print(\"MP2 Corr: {:16.10f}\".format(ref_mp2_corr))\n",
    "print(\"MP3 Corr: {:16.10f}\".format(ref_mp3_corr))\n",
    "print(\"MP4 S   : {:16.10f}\".format(ref_mp4_S))\n",
    "print(\"MP4 D   : {:16.10f}\".format(ref_mp4_D))\n",
    "print(\"MP4 T   : {:16.10f}\".format(ref_mp4_T))\n",
    "print(\"MP4 Q   : {:16.10f}\".format(ref_mp4_Q))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RMP2 相关能量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里基本参照了 Helgaker 书的叙述思路，因此可能与我曾经写过的文档记号和变量名有微妙的区别。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| 变量名 | 公式表达式 | 意义 | 程序的角标顺序 | 出处 |\n",
    "|:---:|:---:|:---:|:---:|:---:|\n",
    "| `E_iajb` | $\\varepsilon_{ij}^{ab}$ | 轨道能差 | $(i, a, j, b)$ | eq (14.4.18) |\n",
    "| `g_mo` | $g_{pqrs}$ | 分子轨道 ERI 积分 | $(p, q, r, s)$ | - |\n",
    "| `t_iajb` | ${t_{ij}^{ab}}^{(1)}$ | 一阶激发系数 | $(i, a, j, b)$ | eq (14.4.41) |\n",
    "| `L_mo` | $L_{pqrs}$ | 定义量 | $(p, q, r, s)$ | eq (13.7.15) |\n",
    "| `T_iajb` | ${\\bar t_{ij}^{ab}}^{(1)}$ | 定义量 | $(i, a, j, b)$ | eq (14.4.42) |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上述定义的一些实际表达式为\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\varepsilon_{ij}^{ab} &= - \\varepsilon_i + \\varepsilon_a - \\varepsilon_j + \\varepsilon_b \\\\\n",
    "g_{pqrs} &= C_{\\mu p} C_{\\nu q} (\\mu \\nu | \\kappa \\lambda) C_{\\kappa r} C_{\\lambda s} \\\\\n",
    "{t_{ij}^{ab}}^{(1)} &= \\frac{g_{iajb}}{\\varepsilon_{ij}^{ab}} \\\\\n",
    "L_{pqrs} &= 2 g_{pqrs} - g_{psrq} \\\\\n",
    "{\\bar t_{ij}^{ab}}^{(1)} &= 2 \\frac{L_{iajb}}{\\varepsilon_{ij}^{ab}} = 4 {t_{ij}^{ab}}^{(1)} - 2 {t_{ij}^{ba}}^{(1)}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "eri_ao = mol.intor(\"int2e\")\n",
    "eri_mo = np.einsum(\"uvkl, up, vq, kr, ls -> pqrs\", eri_ao, C, C, C, C)\n",
    "E_iajb = - e[so, None, None, None] + e[None, sv, None, None] - e[None, None, so, None] + e[None, None, None, sv]\n",
    "g_mo = eri_mo\n",
    "t_iajb = - g_mo[so, sv, so, sv] / E_iajb\n",
    "L_mo = 2 * g_mo - g_mo.swapaxes(-1, -3)\n",
    "T_iajb = 4 * t_iajb - 2 * t_iajb.swapaxes(-1, -3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这些定义下，我们可以相当容易地写出 RMP2 的能量 (Helgaker, eq (14.4.55))\n",
    "\n",
    "$$\n",
    "E_\\mathrm{RMP2, corr} = {t_{ij}^{ab}}^{(1)} L_{iajb}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.26901177170795454"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "energy_mp2_corr = (t_iajb * L_mo[so, sv, so, sv]).sum()\n",
    "energy_mp2_corr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RMP3 相关能"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "RMP3 能量可以根据 Helgaker, eq (14.4.61-62) 构建：\n",
    "\n",
    "$$\n",
    "X_{ij}^{ab} =\n",
    "\\frac{1}{2} {t_{ij}^{cd}}^{(1)} g_{acbd} + \\frac{1}{2} {t_{kl}^{ab}}^{(1)} g_{kilj} + {t_{ik}^{ac}}^{(1)} L_{bjkc} - {t_{kj}^{ac}}^{(1)} g_{bcki} - {t_{ki}^{ac}}^{(1)} g_{bjkc}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "mp3_X_iajb = (\n",
    "    + 0.5 * np.einsum(\"icjd, acbd -> iajb\", t_iajb, g_mo[sv, sv, sv, sv])\n",
    "    + 0.5 * np.einsum(\"kalb, kilj -> iajb\", t_iajb, g_mo[so, so, so, so])\n",
    "    + np.einsum(\"iakc, bjkc -> iajb\", t_iajb, L_mo[sv, so, so, sv])\n",
    "    - np.einsum(\"kajc, bcki -> iajb\", t_iajb, g_mo[sv, sv, so, so])\n",
    "    - np.einsum(\"kaic, bjkc -> iajb\", t_iajb, g_mo[sv, so, so, sv])\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "关于 Helgaker 书上的 eq (14.4.59)，可能是由于我们的实现与之实现方式不同，因此 Helgaker 书中的 ${\\tilde t_{ij}^{ab}}^{(1)}$ 可以当作 ${\\bar t_{ij}^{ab}}^{(1)}$ 使用。因此，\n",
    "\n",
    "$$\n",
    "E_\\mathrm{RMP3, corr} = {\\bar t_{ij}^{ab}}^{(1)} X_{ij}^{ab}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.007451743819516383"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "energy_mp3_corr = (T_iajb * mp3_X_iajb).sum()\n",
    "energy_mp3_corr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RMP4(SDQ) 相关能"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在 RMP4 相关能计算过程中，由于其中涉及到的三激发项从实现上较为复杂，因此我们将会对其分开讨论。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RMP4(S) 相关能"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "定义以下变量：\n",
    "\n",
    "| 变量名 | 公式表达式 | 意义 | 程序的角标顺序 | 出处 |\n",
    "|:---:|:---:|:---:|:---:|:---:|\n",
    "| `E_ia` | $\\varepsilon_{i}^{a}$ | 轨道能差 | $(i, a)$ | eq (14.4.18) |\n",
    "| `t_2_ia` | ${t_{i}^{a}}^{(2)}$ | 二阶激发系数 | $(i, a)$ | eq (14.4.50) |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\varepsilon_i^a$ 的定义是显然的；${t_{i}^{a}}^{(2)}$ 的定义为\n",
    "\n",
    "$$\n",
    "{t_{i}^{a}}^{(2)} = \\frac{1}{\\varepsilon_i^a} \\left[ {t_{kl}^{ad}}^{(1)} L_{kild} - {t_{ki}^{cd}}^{(1)} L_{adkc} \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "定义以下变量：\n",
    "\n",
    "| 变量名 | 公式表达式 | 意义 | 程序的角标顺序 | 出处 |\n",
    "|:---:|:---:|:---:|:---:|:---:|\n",
    "| `E_ia` | $\\varepsilon_{i}^{a}$ | 轨道能差 | $(i, a)$ | eq (14.4.18) |\n",
    "| `t_2_ia` | ${t_{i}^{a}}^{(2)}$ | 二阶激发系数 | $(i, a)$ | eq (14.4.50) |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "E_ia = - e[so, None] + e[None, sv]\n",
    "t_2_ia = (\n",
    "    + np.einsum(\"kald, kild -> ia\", t_iajb, L_mo[so, so, so, sv])\n",
    "    - np.einsum(\"kcid, adkc -> ia\", t_iajb, L_mo[sv, sv, so, sv])\n",
    ")\n",
    "t_2_ia /= E_ia"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由此，根据 Helgaker eq (14.4.79) 与 eq (14.4.83)，有\n",
    "\n",
    "$$\n",
    "S_{ij}^{ab} = {t_{j}^{c}}^{(2)} g_{aibc} - {t_{k}^{b}}^{(2)} g_{aikj}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mp4_S_iajb = (\n",
    "    + np.einsum(\"jc, aibc -> iajb\", t_2_ia, g_mo[sv, so, sv, sv])\n",
    "    - np.einsum(\"kb, aikj -> iajb\", t_2_ia, g_mo[sv, so, so, so])\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以及\n",
    "\n",
    "$$\n",
    "E_\\mathrm{RMP4, S} = {\\bar t_{ij}^{ab}}^{(1)} S_{ij}^{ab}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.004382160888639593"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "energy_mp4_S = (T_iajb * mp4_S_iajb).sum()\n",
    "energy_mp4_S"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RMP4(D) 相关能"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "定义以下变量：\n",
    "\n",
    "| 变量名 | 公式表达式 | 意义 | 程序的角标顺序 | 出处 |\n",
    "|:---:|:---:|:---:|:---:|:---:|\n",
    "| `t_2_iajb` | ${t_{ij}^{ab}}^{(2)}$ | 二阶激发系数 | $(i, a, j, b)$ | eq (14.4.51) |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "${t_{ij}^{ab}}^{(2)}$ 的定义为\n",
    "\n",
    "$$\n",
    "{t_{ij}^{ab}}^{(2)} = \\frac{1}{\\varepsilon_{ij}^{ab}} \\left[ - {t_{ij}^{cd}}^{(1)} g_{acbd} - {t_{kl}^{ab}}^{(1)} g_{kilj} - \\hat P_{ij}^{ab} \\left( {t_{ik}^{ac}}^{(1)} L_{bjkc} - {t_{kj}^{ac}}^{(1)} g_{bcki} - {t_{ki}^{ac}}^{(1)} g_{bjkc} \\right) \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "注意到这里的算符 $\\hat P_{ij}^{ab}$ 是对称化算符，其出处是 Helgaker eq (13.7.13)：\n",
    "\n",
    "$$\n",
    "\\hat P_{ij}^{ab} = A_{ij}^{ab} + A_{ji}^{ba}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为此，我们将 ${t_{ij}^{ab}}^{(2)}$ 的后一个关于 $(i, a, j, b)$ 的张量 $\\left( {t_{ik}^{ac}}^{(1)} L_{bjkc} - {t_{kj}^{ac}}^{(1)} g_{bcki} - {t_{ki}^{ac}}^{(1)} g_{bjkc} \\right)$ 先使用变量 `tmp_symm` 储存，随后使用 `np.transpose` 转置来执行对称化。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp_symm = (\n",
    "    + np.einsum(\"iakc, bjkc -> iajb\", t_iajb, L_mo[sv, so, so, sv])\n",
    "    - np.einsum(\"kajc, bcki -> iajb\", t_iajb, g_mo[sv, sv, so, so])\n",
    "    - np.einsum(\"kaic, bjkc -> iajb\", t_iajb, g_mo[sv, so, so, sv])\n",
    ")\n",
    "tmp_symm += tmp_symm.transpose((2, 3, 0, 1))\n",
    "t_2_iajb = (\n",
    "    - np.einsum(\"icjd, acbd -> iajb\", t_iajb, g_mo[sv, sv, sv, sv])\n",
    "    - np.einsum(\"kalb, kilj -> iajb\", t_iajb, g_mo[so, so, so, so])\n",
    "    - tmp_symm\n",
    ")\n",
    "t_2_iajb /= E_iajb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由此，根据 Helgaker eq (14.4.80) 与 eq (14.4.83)，有\n",
    "\n",
    "$$\n",
    "D_{ij}^{ab} = \\frac{1}{2} {t_{ij}^{cd}}^{(2)} g_{acbd} + \\frac{1}{2} {t_{kl}^{ab}}^{(2)} g_{kilj} + {t_{ik}^{ac}}^{(2)} L_{bjkc} - {t_{kj}^{ac}}^{(2)} g_{bcki} - {t_{ki}^{ac}}^{(2)} g_{bjkc}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "mp4_D_iajb = (\n",
    "    + 0.5 * np.einsum(\"icjd, acbd -> iajb\", t_2_iajb, g_mo[sv, sv, sv, sv])\n",
    "    + 0.5 * np.einsum(\"kalb, kilj -> iajb\", t_2_iajb, g_mo[so, so, so, so])\n",
    "    + 1.0 * np.einsum(\"iakc, bjkc -> iajb\", t_2_iajb, L_mo[sv, so, so, sv])\n",
    "    - 1.0 * np.einsum(\"kajc, bcki -> iajb\", t_2_iajb, g_mo[sv, sv, so, so])\n",
    "    - 1.0 * np.einsum(\"kaic, bjkc -> iajb\", t_2_iajb, g_mo[sv, so, so, sv])\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以及\n",
    "\n",
    "$$\n",
    "E_\\mathrm{RMP4, D} = {\\bar t_{ij}^{ab}}^{(1)} D_{ij}^{ab}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.008819028049601586"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "energy_mp4_D = (T_iajb * mp4_D_iajb).sum()\n",
    "energy_mp4_D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RMP4(Q) 相关能"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在 Helgaker 的书中，RMP4(Q) 的相关能没有再作额外定义。根据 Helgaker eq (14.4.82) 与 eq (14.4.83)，有"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "Q_{ij}^{ab} &=\n",
    "  \\frac{1}{2} {t_{kl}^{ab}}^{(1)} {t_{ij}^{cd}}^{(1)} g_{kcld}\n",
    "+             {t_{ik}^{ac}}^{(1)} {t_{jl}^{bd}}^{(1)} L_{kcld}\n",
    "-             {t_{ik}^{ac}}^{(1)} {t_{lj}^{bd}}^{(1)} L_{kcld} \\\\ & \\quad\n",
    "+ \\frac{1}{2} {t_{ki}^{ac}}^{(1)} {t_{lj}^{bd}}^{(1)} g_{kcld}\n",
    "+ \\frac{1}{2} {t_{kj}^{ad}}^{(1)} {t_{li}^{bc}}^{(1)} g_{kcld} \\\\ & \\quad\n",
    "-             {t_{ik}^{ab}}^{(1)} {t_{lj}^{cd}}^{(1)} L_{lckd}\n",
    "-             {t_{ij}^{ac}}^{(1)} {t_{kl}^{bd}}^{(1)} L_{kcld}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "mp4_Q_iajb = (\n",
    "    + 0.5 * np.einsum(\"kalb, icjd, kcld -> iajb\", t_iajb, t_iajb, g_mo[so, sv, so, sv])\n",
    "    + 1.0 * np.einsum(\"iakc, jbld, kcld -> iajb\", t_iajb, t_iajb, L_mo[so, sv, so, sv])\n",
    "    - 1.0 * np.einsum(\"iakc, lbjd, kcld -> iajb\", t_iajb, t_iajb, L_mo[so, sv, so, sv])\n",
    "    + 0.5 * np.einsum(\"kaic, lbjd, kcld -> iajb\", t_iajb, t_iajb, g_mo[so, sv, so, sv])\n",
    "    + 0.5 * np.einsum(\"kajd, lbic, kcld -> iajb\", t_iajb, t_iajb, g_mo[so, sv, so, sv])\n",
    "    - 1.0 * np.einsum(\"iakb, lcjd, lckd -> iajb\", t_iajb, t_iajb, L_mo[so, sv, so, sv])\n",
    "    - 1.0 * np.einsum(\"iajc, kbld, kcld -> iajb\", t_iajb, t_iajb, L_mo[so, sv, so, sv])\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以及\n",
    "\n",
    "$$\n",
    "E_\\mathrm{RMP4, Q} = {\\bar t_{ij}^{ab}}^{(1)} Q_{ij}^{ab}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.001135577514293328"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "energy_mp4_Q = (T_iajb * mp4_Q_iajb).sum()\n",
    "energy_mp4_Q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由此，我们可以计算 RMP4(SDQ) 所贡献的相关能："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.01206561142394785"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "energy_mp4SDQ_corr = energy_mp4_S + energy_mp4_D + energy_mp4_Q\n",
    "energy_mp4SDQ_corr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RMP4(T) 相关能"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "RMP4(T) 从实现上来讲，最简单的做法需要消耗 $O^3 V^3$ 的内存。我们先从简单的实现入手；随后会简要给出一个中间张量的内存消耗是 $O^2 V^2$ 的算法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RMP4(T) 相关能：简单实现"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "定义以下变量：\n",
    "\n",
    "| 变量名 | 公式表达式 | 意义 | 程序的角标顺序 | 出处 |\n",
    "|:---:|:---:|:---:|:---:|:---:|\n",
    "| `E_iajbkc` | $\\varepsilon_{ijk}^{abc}$ | 轨道能差 | $(i, a, j, b, k, c)$ | eq (14.4.18) |\n",
    "| `t_2_iajbkc` | ${t_{ijk}^{abc}}^{(2)}$ | 二阶激发系数 | $(i, a, j, b, k, c)$ | eq (14.4.52) |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\varepsilon_{ijk}^{abc}$ 的定义是显然的；${t_{ijk}^{abc}}^{(2)}$ 的定义为\n",
    "\n",
    "$$\n",
    "{t_{ijk}^{abc}}^{(2)} = - \\frac{1}{\\varepsilon_{ijk}^{abc}} \\hat P_{ijk}^{abc} \\left( {t_{ij}^{ad}}^{(1)} g_{ckbd} - {t_{il}^{ab}}^{(1)} g_{cklj} \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里的算符 $\\hat P_{ijk}^{abc}$ 仍然是对称算符，但其形式更复杂 (Helgaker, eq (13.7.14))：\n",
    "\n",
    "$$\n",
    "\\hat P_{ijk}^{abc} A_{ijk}^{abc} = A_{ijk}^{abc} + A_{ikj}^{acb} + A_{jik}^{bac} + A_{jki}^{bca} + A_{kij}^{cab} + A_{kji}^{cba}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "E_iajbkc = (\n",
    "    - e[so, None, None, None, None, None]\n",
    "    + e[None, sv, None, None, None, None]\n",
    "    - e[None, None, so, None, None, None]\n",
    "    + e[None, None, None, sv, None, None]\n",
    "    - e[None, None, None, None, so, None]\n",
    "    + e[None, None, None, None, None, sv]\n",
    ")\n",
    "t_2_iajbkc = (\n",
    "    + np.einsum(\"iajd, ckbd -> iajbkc\", t_iajb, g_mo[sv, so, sv, sv])\n",
    "    - np.einsum(\"ialb, cklj -> iajbkc\", t_iajb, g_mo[sv, so, so, so])\n",
    ")\n",
    "t_2_iajbkc = (\n",
    "    + t_2_iajbkc.transpose((0, 1, 2, 3, 4, 5))\n",
    "    + t_2_iajbkc.transpose((0, 1, 4, 5, 2, 3))\n",
    "    + t_2_iajbkc.transpose((2, 3, 0, 1, 4, 5))\n",
    "    + t_2_iajbkc.transpose((2, 3, 4, 5, 0, 1))\n",
    "    + t_2_iajbkc.transpose((4, 5, 0, 1, 2, 3))\n",
    "    + t_2_iajbkc.transpose((4, 5, 2, 3, 0, 1))\n",
    ")\n",
    "t_2_iajbkc = - t_2_iajbkc / E_iajbkc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由此，根据 Helgaker eq (14.4.81) 与 eq (14.4.83)，有\n",
    "\n",
    "$$\n",
    "T_{ij}^{ab} = \n",
    "  {t_{ijk}^{acd}}^{(2)} L_{bckd}\n",
    "- {t_{kji}^{acd}}^{(2)} g_{kdbc}\n",
    "- {t_{ikl}^{abc}}^{(2)} L_{kjlc}\n",
    "+ {t_{lki}^{abc}}^{(2)} g_{kjlc}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "mp4_T_iajb = (\n",
    "    + np.einsum(\"iajckd, bckd -> iajb\", t_2_iajbkc, L_mo[sv, sv, so, sv])\n",
    "    - np.einsum(\"kajcid, kdbc -> iajb\", t_2_iajbkc, g_mo[so, sv, sv, sv])\n",
    "    - np.einsum(\"iakblc, kjlc -> iajb\", t_2_iajbkc, L_mo[so, so, so, sv])\n",
    "    + np.einsum(\"lakbic, kjlc -> iajb\", t_2_iajbkc, g_mo[so, so, so, sv])\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以及\n",
    "\n",
    "$$\n",
    "E_\\mathrm{RMP4, T} = {\\bar t_{ij}^{ab}}^{(1)} T_{ij}^{ab}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.006712603570763939"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "energy_mp4_T = (T_iajb * mp4_T_iajb).sum()\n",
    "energy_mp4_T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由此，我们可以给出完整的 RMP4 相关能矫正了："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.018778214994711787"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "energy_mp4_corr = energy_mp4_S + energy_mp4_D + energy_mp4_T + energy_mp4_Q\n",
    "energy_mp4_corr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RMP4(T) 相关能：中间张量 $O^2 V^2$ 内存实现"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里我们通过限制 ${t_{ijk}^{abc}}^{(2)}$ 中的 $b$ 与 $k$ 并作求和的方式，将中间矩阵的储存大小降为 $O^2 V^2$ 并给出 RMP4(T) 的能量。\n",
    "\n",
    "这个严格来说还是至少消耗了 $O^1 V^3$ 的内存，因为在计算过程中用到了张量 $g_{aicd}$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "首先指出，像上述代码中使用 `np.transpose` 的方式来处理 $\\hat P_{ijk}^{abc}$ 在这里可能不适用；因此需要手输所有的中间张量的缩并过程。下述函数可以用来辅助我们进行带有指标转换的张量缩并。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://stackoverflow.com/a/11122744/9647779\n",
    "def substitute_string(string, rule):\n",
    "    str_lst = list(string)\n",
    "    rule1, rule2 = rule.replace(\" \", \"\").split(\"->\")\n",
    "    idx_list = [[i.start() for i in re.finditer(c, string)] for c in rule1]\n",
    "    for idx, pos_list in enumerate(idx_list):\n",
    "        for pos in pos_list:\n",
    "            str_lst[pos] = rule2[idx]\n",
    "    return \"\".join(str_lst)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "譬如我们现在需要将缩并的目标 ${t_{ij}^{ad}}^{(1)} g_{ckbd}$ 的 $(i, a, j, b, k, c)$ 转换为 $(j, b, k, c, i, a)$，那么我们执行下述代码就可以生成转换后的缩并字符串了："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'jbkd, aicd'"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "substitute_string(\"iajd, ckbd\", \"iajbkc -> jbkcia\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "后文中出现的代码就是受这个小函数的帮助而写成的。\n",
    "\n",
    "事实上，下述程序几乎与上面的程序等价，只是每次处理 ${t_{ijk}^{abc}}^{(2)}$ 张量时，总是先固定 $b, k$ 指标，而对其它的张量进行正常的运算，求出固定了 $b, k$ 的 RMP4(T) 相关能；随后我们再对 $b, k$ 指标进行循环，把所有 $b, k$ 相关能贡献总和起来，得到最终的 RMP4(T) 相关能。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.006712603570763938"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "energy_mp4_T_by_ckiter = 0\n",
    "for k in range(nocc):\n",
    "    for b in range(nvir):\n",
    "        bn = b + nocc\n",
    "        t_2_bk_iajb = - (\n",
    "            + np.einsum(\"iajd, cd -> iajc\", t_iajb[:, :, :, :], g_mo[sv,  k, bn, sv])\n",
    "            - np.einsum(\"ial, clj -> iajc\", t_iajb[:, :, :, b], g_mo[sv,  k, so, so])\n",
    "            + np.einsum(\"iad, jcd -> iajc\", t_iajb[:, :, k, :], g_mo[bn, so, sv, sv])\n",
    "            - np.einsum(\"ialc, jl -> iajc\", t_iajb[:, :, :, :], g_mo[bn, so, so,  k])\n",
    "            + np.einsum(\"jid, cad -> iajc\", t_iajb[:, b, :, :], g_mo[sv,  k, sv, sv])\n",
    "            - np.einsum(\"jla, cli -> iajc\", t_iajb[:, b, :, :], g_mo[sv,  k, so, so])\n",
    "            + np.einsum(\"jd, aicd -> iajc\", t_iajb[:, b, k, :], g_mo[sv, so, sv, sv])\n",
    "            - np.einsum(\"jlc, ail -> iajc\", t_iajb[:, b, :, :], g_mo[sv, so, so,  k])\n",
    "            + np.einsum(\"cid, jad -> iajc\", t_iajb[k, :, :, :], g_mo[bn, so, sv, sv])\n",
    "            - np.einsum(\"cla, jli -> iajc\", t_iajb[k, :, :, :], g_mo[bn, so, so, so])\n",
    "            + np.einsum(\"cjd, aid -> iajc\", t_iajb[k, :, :, :], g_mo[sv, so, bn, sv])\n",
    "            - np.einsum(\"cl, ailj -> iajc\", t_iajb[k, :, :, b], g_mo[sv, so, so, so])\n",
    "        )\n",
    "        t_2_bk_iajb /= (\n",
    "            - e[so, None, None, None]\n",
    "            + e[None, sv, None, None]\n",
    "            - e[None, None, so, None]\n",
    "            + e[None, None, None, sv]\n",
    "            - e[k] + e[bn]\n",
    "        )\n",
    "        mp4_T_bk_iajb = (\n",
    "            + np.einsum(\"iajd, bd -> iajb\", t_2_bk_iajb, L_mo[sv, bn,  k, sv])\n",
    "            - np.einsum(\"idja, db -> iajb\", t_2_bk_iajb, g_mo[ k, sv, sv, bn])\n",
    "            - np.einsum(\"ialb, jl -> iajb\", t_2_bk_iajb, L_mo[ k, so, so, bn])\n",
    "            + np.einsum(\"laib, jl -> iajb\", t_2_bk_iajb, g_mo[ k, so, so, bn])\n",
    "        )\n",
    "        energy_mp4_T_by_ckiter += (T_iajb * mp4_T_bk_iajb).sum()\n",
    "energy_mp4_T_by_ckiter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[^Helgaker-Jorgensen.Wiley.2013]: Helgaker, T.; Olsen, J.; Jorgensen, P. *Molecular Electronic-Structure Theory*; Wiley-Blackwell, 2013-02-01.\n",
    "\n",
    "[^Szabo-Ostlund.Dover.1996]: Szabo, A.; Ostlund, N. S. *Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory (Dover Books on Chemistry)*; Dover Publications, 1996.\n",
    "\n",
    "[^psi4numpy-mp3]: <https://github.com/psi4/psi4numpy/tree/master/Moller-Plesset>"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
